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Abstract 

The relaxation of polarized spins in a porous medium has been utilized as a probe 
of its structure. We note that the governing diffusion problem has a close parallel to 
that of a particle in a box, an elementary Quantum mechanics toy model. Following 
the spirits of "free electron" model, we use generic properties of the eigen spectrum 
to understand features common to a wide variety of pore geometry, consistent with 
large scale numerical simulations and experimental data. 
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INTRODUCTION 

Relaxation of polarized spins by paramagnetic impurities embedded in pore matrix provides 
a valuable probe into its structure, especially when elaborate techniques such as MRI imag- 
ing are not easily applicable. Exploration for hydrocarbon inside a bore hole is one such 
case where it is critical to extract limited, but pertinent information in a timely manner 
from measurements often done with a hand tied. Over the past years, a lot of progress has 
been made both in development of magnetic resonance tools 1 as well as theoretical aspects 2 
geared specifically toward extracting surface-to-volume ratio from the observed magnetiza- 
tion evolution. There are several more sophisticated MR techniques now being routinely 
applied to rocks in the laboratories and thus potentially much more detailed information on 
the structure of the pore space is available. 3 In addition to the pore structure which con- 
trols practically important properties such as ite permeability, one is often interested in the 
details of pore filling fluids or gas such as how is is made up out of water and oil, diffusion 
constants of its constituents, and wetting condition. As our scope begins to embrace this 
degree of complexity, it becomes essential to bring some concepts up to date since sorting 
out the complexity of the pore matrix from that of its content is a minimum requirement in 
such endeavor. 

Although we are specifically motivated by MR applied to porous media, the underlying 
diffusion phenemenon under complex boundary conditions occurs in a variety of real life 
situations. To illustrate a useful point, let us consider an extreme case: zombies wandering 
on top of a shopping mall whose fence has occasional openings at the edge. Depending on 
the size of the space, the speed of a zombie, and the frequency of the openings, the depletion 
of its population occurs either in a homogeneous way or with more interesting profile in 
space. The conventional approach sans imaging relied on deducing these parameters from 
observing how the overall population decays. 

QUANTUM MECHANICS 

The analysis is simple if we take the continuum approach, treating the population as a 
scalar variable, subject to the classic diffusion equation. 4 ' 5 The diffusive current jo arises 
from inhomogeneous density distribution m via ju = —DVin. Conservation of local flux 
V • jn = ~^i m then leads to DX7 2 m = j^m with the boundary condition h ■ jd = pm on the 
wall. Treating it as a boundary value problem, the problem is solved in terms of eigenmode 
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expansion, as was studied in detail by Brownstein and Taar for simple geometries years 
earlier. In this work, we employ an analogy to Quantum mechanics by observing that 
D — > y~ and t — > it turns the equation for m into the Schrodinger eqaution for a particle 
confined to a leaky container with Ti = — ^V 2 . m is identified with the wave function ^ 
of the QM particle in an imaginary time, but we note that the m is directly observed in 
our classical system, rather than |\P| 2 . The Hamiltonian operator is simply that of kinetic 
energy of the particle and the boundary condition is of mixed type, [n ■ DV + p]m = 0. The 
existence, orthonormality of eigenmodes are guarranteed with satisfying H^ n = E n ^ n . 
By multiplying both sides with \&* and integrating over the pore volume, with the help of 
the boundary condition, one can show 

{^J dh-\V^ n \ 2 + phj\ds\\m n \ 2 ) =E n , (0.1) 

which shows that E n > with equality allowed only for p = with a uniform mode = 
The energy (relaxation rate) is broken down into two components: the kinetic energy part 
and the potential energy part which is localized in the surface region. For each eigenmode, 
we have a length scale £ n = $-^^r ,o with which we have 

En= «_j dW(l + ^ tn) . ,0.2) 

Since £ n oc 4 A- decreases with increasing n, this clearly shows how the high (excited) modes 
become progressively independent of p. 

From now on, we introduce the bra-ket notation, \l/ n (r) =< v\n > for compactness. The 
eigenmodes form a complete set with \n > J2^=o < n\ = 1 and the time evolution of an 
arbitrary state > is given by 

\^ t >=^2\n > e'^ < n\^i > . (0.3) 

n 

This is all very formal. The hard part lies in the specific spectrum for a given boundary 
shape and the strength of the sink. The distribution of eigenvalues encodes much informa- 
tion of the boundary shape, i.e. pore structure. 6 It is possible that there exists significant 
degeneracy, since different shapes may lead to an identical spectrum. Despite this, we an- 
ticipate that we can relate some characteristics of the eigenspectrum to a given shape of 
the pore wall. Borrowing the concepts from random matrix theory and Quantum chaos, 
we may study the correlation among the levels via evaluating V(e nm ) where e = E n — E m . 
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This function generally displays a universal characterisitc for large E n 's, while at small val- 
ues (correlation among slow diffusion modes) unique fingerprint of the boundary structure 
manifests itself. Note that p acts as a knob taking the system from a Neunmann (p — > 0) 
to a Dirchlet (p — ► oo) boundary condition. This variation in sink strength largely affect 
the lowest eigen mode in a non-trivial manner while for higher levels, amounts only to a 
uniform phase shift by ~ ~. The effect of a large scale pore shape on the spectrum can be 
most easily appreciated in the model system of particle in a rectangle. There the spectrum 
follows ( in the p — > limit for simplicity): 

^ = ^ = h f[Cr? + C^) 2 + Crn (o.4) 

fl 1 n ZTYl LJ X Ajy Lj z 

One can easilty obtain the density of states from the surface area of an ellisoid in the k 
space as in the free electron model. Some representative spectra are shown in Figure |2| in 
which we varied the aspect ratio of a rectangle while keeping its volume fixed at L 3 = 1 
and p = 0. Qualitatively, the shape of the bouding boxes represents a 1D-, 2D- and 3D 
space respectively and the corresponding spectrum displays a characteristic distribution. As 
mentioned earlier, what is probably more meaningful in respect of the boundary shape is 
the auto-correlation function of the level positions. Practically, this requires that a given 
measurement picks out the full spectrum with extremely high resolution. Unfortunately, 
this step is highly compromised by the inherent limitations of inverse Laplace transform 
with which one transforms the time-domain data into the spectrum. In addition to this, 
the preparation and detection states usually have high degree of symmetry, thus picking out 
only a subset of the full spectrum. 

APPLICATION 

Usually, the final state is detected by measuring the net magnetization. This is equivalent 
to projecting the final state onto a uniform state. If one also prepared a uniform initial state, 
l^j >= |^o >> then we obtain 

M(t) =< ^o|^ >t=Y, e ~ Ett \ < n \^o > I 2 - (0.5) 



For small p, we know from Equation |0J] that the lowest eigenmode becomes close to the 
uniform mode \I/o ~ with E ~ pS/V following directly from the equation. Furthermore, 
from orthogonality condition, we expect < n|0 >= 5 n o and therefore, A4(t) ~ e~~^ 1 . This 
immediately shows the limitations of detecting uniform magnetization as a probe for multiple 



length scales in a given pore. Even for determination of the average surface-to-volume ratio, 
y from the observed E ~ py is problematic unless the value of p is independently and 
accurately known since the often dominant, lowest eigenmode is extremely sensitive to p 
and the all higher modes which remain robust as p is varied, unfortunately do not show up 
in measurements. To illustrate this quantitatively, we show in Figure ^ partial spectrum 
for a spherical pore of radius a as we vary the control parameter pa/D from 0.05 to 5. 
The general eigenmodes for a sphere is given by ^k,i,m — jiihrfY™ with I = 0, 1, . . . and 
m = —I, ... 0, ... /. For given I and m, there is a manifold of eigenmodes with k determied 
from 

fca ^-i(^)-^ + +l^( fcfl ) + ^ji(ka) = 0. (0.6) 

Usually, only the subset with I = m = is considered on symmetrical grounds with isotropic 
initial and detected states >> but for more general cases (See discussion below), full 
angular variation should be retained even with the isotropic boundary condition. The most 
salient feature in the figure is the extreme sensitivity of the lowest eigenmode on p. The 
higher modes, including the lowest modes in each I ^ manifolds, are less sensitive and 
their positions give direct information on the size of the sphere a. Another important 
characteristics of the higher modes is that they generally lead to a universal density of 
states, a la free electron model, which reflects the fact that the most pore walls look very 
similar on microscopic length scales, and the fast decay modes are insensitive to the large 
scale variation of the pore wall. 

There are ways to get around the barrier to the excited states. One way is to create and 
detect a localized state. This extreme case begins with a localized polarization < r|\E^ >= 
S(r — ro), and detecting again at the same spot, < The observed signal will follow: 

M loc (r ) =J2e- t/T -\ < n\r > \ 2 (0.7) 

n 

with the spectral weight | < ro|n > | 2 for each mode democratically distributed over a 
broad range. The two cases considered so far are two extreme limits where the prepared and 
the detected states are oriented in the Hilbert space with some overlap with the (uniform) 
ground state of the given boundary value problem, as schematically depicted in Figure |j. 
The most general situation may be described in terms of a unitary rotation of the uniform 
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state by the application of a local phase encoding-decoding operation, e ia ^\ 

M a = J2e~ t/Tn \ < n|e M(r) |^ > ? ■ (0.8) 

n 

This operation can indeed be realized in the laboratory via application of 90- and 180-degree 
rf pulses, a common MR techque. The critical ingredient for this to work is the presence of 
internal field inhomogeneity on sub-pore length scales. In the presence of a uniform applied 
field along the z-axis, the small magnetic susceptibility contrast between the matrix and its 
filling fluid induces a local field proportional to A%. Therefore, one may rotate the spins 
onto the plane transverse to the polarizing field, and let them precess under the local field 
h(r) for a given time r, accumulating a phase 0(r) = 7/1(1:) t. 

The spins are then rotated back, orthogonal to the transverse plane and allowed to relax 
over time t before one apply the reversing process. This winding-unwinding procedure would 
work perfectly if the molecules do not diffuse. However, the wound (or encoded) state now 
evolves with \n > e~ l l Tn < n\ with significantly enhanced overlap with all excited states 
\n > with n > 0. One can immediately see a potential benefit if the modulating field h(r) 
closely follows the underlying geometry: The rotation of > will occur in such a way that 
the maximal overalp with the class of eigenmodes with length scales matching the variation 
of the boundary One can qualitatively understand this by noting that the local field can be 
described as arising from an inhomogeneous distribution of surface dipoles whose strength 
is proportional to h-n. Therefore, the local field tends to get enhanced in the vicinity of the 
pore wall, and is sensitive to its local orientation. The encoded pattern created by flashing of 
the local field (in practice, the spins are momentarily rotated to feel its presence) therefore 
reflects the contour of the pore wall in a convoluted manner. Notwithstanding some detailed 
issues, this is quite sufficient in bringing out the full eigenspectrum into light. 

Finally, we note that the widely studied D(t) oc \r(t) — r(0) | 2 is analogous to the polar- 
izability J2 n < ^ f\r\n > e~ tn < n\r\^i > which for a rectangle leads to ^{L 2 X + L 2 y + L 2 Z ) + 

, DniT 2 An 2 
-t " . 2 j_2 

^2 J2 n >o Si e L% n?(.~ !)""• This, along with generic properties of the excited states, leads 
to an analytic for D(t) which covers a wide range of t turning it into an effeftive probe into 
the pore shape as will be reported elsewhere. 

To be concrete, we are carrying out large scale numerical simulations of spin diffusion in 
various pore structures such as random glass bead pack and real rocks obtained with 3D 
microtomographic technique. 8 The complex 3D pores are represented as coupled cluster of 
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micron length scale cubes in which a large number of random walkers (typically 10 4 — 10 5 in 
numbers) are deployed to simulate a variety of MR measurements. For A4 a , we obtain the 
internal field by numerically evaluating dipole contributions from all surrounding chambers. 
Looking at the results, an example of which is shown in Figure |5], one can indeed verify the 
claim made earlier on the connection between the pore shape and the local field configuration. 

Figure ^] shows a typical simulation results performed on the 3D microtomographic data 
of Fontainebleu sample. The result for M. shows excellent agreement with experimental 
data on the same class of rocks. Furthermore, the simulated M. a excavates the length scales 
hidden in Ai as shown in the inset for | < n|e m ^ r - ) |\[ / o > | 2 - The presence of multiple length 
scales manifests itself in apprearence of the peaks, and their spectral weight depends on the 
strength of a oc AxBt. More detailed numerical studies on the variety of pore structure is 
under progress and will be reported elsewhere. 
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FIG. 1: Evolution of a population of random walkers subject to sinks on the edge. 
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FIG. 2: Spectral "fingerprints" of the boundary shape: For a given volume L 3 (porosity), we com- 
pare three distinct aspect ratio ,^,%) = (11.11, 0.3, 0.3) (top:rod), (1.826,1.826,0.3) (mid- 
dle:slab) and (1,1,1) (bottomxube). T* = L^/I?7r 2 :nominal diffusion time along the longest 
length, p = is used. 
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FIG. 3: Spectral variation for a sphere of radius a. The control parameter pa/D is varied from 
0.05 to 5.0. For each value of this, three groups are shown with angular momentum I = 0, 1,2. 
The arrow marks the lowest eigen mode which shows an extreme sensitivity to p. 
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FIG. 4: Comparison between the probe using uniform initial & final states (left) and that using 
states rotated in the Hilbert sub-space spanned by the eigenmodes conforming to the symmetry 
selection rule of the probe operator. 
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FIG. 5: A typical 3D pore strucure used in numerical simulation. The pore (light gray) is composed 
of about 70000 connected cubes of size 5.72 3 /xm 3 . The inset shows a typical cross section with color 
coding representing the local field variation. 
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FIG. 6: Simulated M(t) and A4 a performed on the Fontainebleu sandstone structure. The dots 
are from NMR experiment by Hiirlimann et al. 7 The inset shows the spectral weight derived by 
performing an inverse Laplace transform. 
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